#ifndef  GMATH_INC
#define  GMATH_INC

/** 
 * @file gmath.h
 * @brief Mathematical constants and functions
 * @author Anshul Nigham
 * @date 2007-07-14
 */

#include <gsl/gsl_math.h>

namespace gslcpp {
	/** 
	 * @brief Class wrapping mathematical constants and functions
	 */
	class gmath
	{
		public:
			/** 
			 * @brief The base of exponentials, \f$e\f$
			 */
			static const double E;
			/** 
			 * @brief The base-2 logarithm of \f$e\f$, \f$\log_2 (e)\f$
			 */
			static const double LOG2E;
			/** 
			 * @brief The base-10 logarithm of \f$e\f$, \f$\log_{10} (e)\f$
			 */
			static const double LOG10E;
			/** 
			 * @brief The square root of 2, \f$\sqrt{2}\f$
			 */
			static const double SQRT2;
			/** 
			 * @brief The square root of one-half, \f$\sqrt{1/2}\f$
			 */
			static const double SQRT1_2;
			/** 
			 * @brief The square root of three, \f$\sqrt{3}\f$
			 */
			static const double SQRT3;
			/** 
			 * @brief The constant pi, \f$\pi\f$
			 */
			static const double PI;
			/** 
			 * @brief Pi divided by 2, \f$\pi / 2\f$
			 */
			static const double PI_2;
			/** 
			 * @brief Pi divided by 4, \f$\pi / 4\f$
			 */
			static const double PI_4;
			/** 
			 * @brief The square root of Pi, \f$\sqrt{\pi}\f$
			 */
			static const double SQRTPI;
			/** 
			 * @brief Two divided by the square root of Pi, \f$2 / \sqrt{\pi}\f$
			 */
			static const double D2_SQRTPI;
			/** 
			 * @brief The reciprocal of Pi, \f$1 / \pi\f$
			 */
			static const double D1_PI;
			/** 
			 * @brief Twice the reciprocal of Pi, \f$2 / \pi\f$
			 */
			static const double D2_PI;
			/** 
			 * @brief The natural logarithm of 10, \f$\ln(10)\f$
			 */
			static const double LN10;
			/** 
			 * @brief The natural logarithm of 2, \f$\ln(2)\f$
			 */
			static const double LN2;
			/** 
			 * @brief The natural logarithm of Pi, \f$\ln(\pi)\f$
			 */
			static const double LNPI;
			/** 
			 * @brief Euler's constant, \f$\gamma\f$
			 */
			static const double EULER;

			/** 
			 * @brief The IEEE representation of positive infinity, \f$+\infty\f$. It is computed from the expression +1.0/0.0.
			 */
			static const double POSINF;
			/** 
			 * @brief The IEEE representation of negative infinity, \f$-\infty\f$. It is computed from the expression -1.0/0.0.
			 */
			static const double NEGINF;
			/** 
			 * @brief The IEEE representation of NaN. It is computed from the expression 0.0/0.0.
			 */
			static const double G_NAN;

			/** 
			 * @brief This function returns TRUE if x is not a number.
			 */
			static bool is_nan (double x);
			/** 
			 * @brief Returns +1 if x is positive infinity, -1 if x is negative infinity, and 0 otherwise.
			 */
			static int is_inf (double x);
			/** 
			 * @brief Returns TRUE if x is a real number, FALSE if it is infinite or not-a-number.
			 */
			static bool is_finite (double x);

			/** 
			 * @brief This function computes the value of \f$\log(1+x)\f$ in a way that is accurate for small x.
			 */
			static double log1p (double x);
			/** 
			 * @brief This function computes the value of \f$\exp(x)-1\f$ in a way that is accurate for small x.
			 */
			static double expm1 (double x);
			/** 
			 * @brief This function computes the value of \f$\sqrt{x^2 + y^2}\f$ in a way that avoids overflow.
			 */
			static double hypot (double x, double y);
			/** 
			 * @brief This function computes the value of \f$\mathrm{arccosh}(x)\f$.
			 */
			static double acosh (double x);
			/** 
			 * @brief This function computes the value of \f$\mathrm{arcsinh}(x)\f$.
			 */
			static double asinh (double x);
			/** 
			 * @brief This function computes the value of \f$\mathrm{arctanh}(x)\f$.
			 */
			static double atanh (double x);
			/** 
			 * @brief This function computes the value of \f$x \times 2^e\f$.
			 */
			static double ldexp (double x, int e);
			/** 
			 * @brief This function splits the number x into its normalized fraction f and exponent e, such that \f$x = f \times 2^e\f$ and \f$0.5 <= f < 1\f$. 
			 */
			static double frexp (double x, int* e);

			/** 
			 * @brief This function computes \f$x^n\f$.
			 */
			static double pow_int (double x, int n);
			/** 
			 * @brief This function computes \f$x^2\f$.
			 */
			static double pow_2 (double x);
			/** 
			 * @brief This function computes \f$x^3\f$.
			 */
			static double pow_3 (double x);
			/** 
			 * @brief This function computes \f$x^4\f$.
			 */
			static double pow_4 (double x);
			/** 
			 * @brief This function computes \f$x^5\f$.
			 */
			static double pow_5 (double x);
			/** 
			 * @brief This function computes \f$x^6\f$.
			 */
			static double pow_6 (double x);
			/** 
			 * @brief This function computes \f$x^7\f$.
			 */
			static double pow_7 (double x);
			/** 
			 * @brief This function computes \f$x^8\f$.
			 */
			static double pow_8 (double x);
			/** 
			 * @brief This function computes \f$x^9\f$.
			 */
			static double pow_9 (double x);

			/** 
			 * @brief This function returns the sign of x. It is defined as ((x) >= 0 ? 1 : -1). Note that with this definition the sign of zero is positive (regardless of its ieee sign bit).
			 */
			static int sign (float x);
			/** 
			 * @brief This function returns the sign of x. It is defined as ((x) >= 0 ? 1 : -1). Note that with this definition the sign of zero is positive (regardless of its ieee sign bit).
			 */
			static int sign (double x);
			/** 
			 * @brief This function returns the sign of x. It is defined as ((x) >= 0 ? 1 : -1). Note that with this definition the sign of zero is positive (regardless of its ieee sign bit).
			 */
			static int sign (int x);
			/** 
			 * @brief This function returns the sign of x. It is defined as ((x) >= 0 ? 1 : -1). Note that with this definition the sign of zero is positive (regardless of its ieee sign bit).
			 */
			static int sign (long x);

			/** 
			 * @brief This function returns TRUE if the integer x is odd.
			 */
			static bool is_odd (int x);
			/** 
			 * @brief This function returns TRUE if the integer x is even.
			 */
			static bool is_even (int x);

			/** 
			 * @brief This function returns the maximum of the numbers x and y.
			 */
			static double max (double x, double y);
			/** 
			 * @brief This function returns the maximum of the numbers x and y.
			 */
			static long double max (long double x, long double y);
			/** 
			 * @brief This function returns the maximum of the numbers x and y.
			 */
			static int max (int x, int y);
			/** 
			 * @brief This function returns the minimum of the numbers x and y.
			 */
			static double min (double x, double y);
			/** 
			 * @brief This function returns the minimum of the numbers x and y.
			 */
			static long double min (long double x, long double y);
			/** 
			 * @brief This function returns the minimum of the numbers x and y.
			 */
			static int min (int x, int y);

			/** 
			 * @brief This function determines whether x and y are approximately equal to a relative accuracy epsilon.
			 *
			 * The relative accuracy is measured using an interval of size \f$2 \delta\f$, where \f$\delta = 2^k \epsilon\f$ and k is the maximum base-2 exponent of x and y as computed by the function frexp.
			 * 
			 * The implementation is based on the package fcmp by T.C. Belding.
			 *
			 * @return Returns 0 if x and y are approximately equal, -1 if x < y, and +1 x > y.
			 */
			static int fcmp (double x, double y, double epsilon);
	};

	const double gmath::E = M_E;
	const double gmath::LOG2E = M_LOG2E;
	const double gmath::LOG10E = M_LOG10E;
	const double gmath::SQRT2 = M_SQRT2;
	const double gmath::SQRT1_2 = M_SQRT1_2;
	const double gmath::SQRT3 = M_SQRT3;
	const double gmath::PI = M_PI;
	const double gmath::PI_2 = M_PI_2;
	const double gmath::PI_4 = M_PI_4;
	const double gmath::SQRTPI = M_SQRTPI;
	const double gmath::D2_SQRTPI = M_2_SQRTPI;
	const double gmath::D1_PI = M_1_PI;
	const double gmath::D2_PI = M_2_PI;
	const double gmath::LN10 = M_LN10;
	const double gmath::LN2 = M_LN2;
	const double gmath::LNPI = M_LNPI;
	const double gmath::EULER = M_EULER;

	const double gmath::POSINF = GSL_POSINF;
	const double gmath::NEGINF = GSL_NEGINF;
	const double gmath::G_NAN = GSL_NAN;

	inline bool gmath::is_nan (double x)
	{
		return gsl_isnan(x);
	}
	inline int gmath::is_inf (double x)
	{
		return gsl_isinf(x);
	}
	inline bool gmath::is_finite ( double x )
	{
		return gsl_finite(x);
	}
	inline double gmath::log1p ( double x )
	{
		return gsl_log1p(x);
	}
	inline double gmath::expm1 ( double x )
	{
		return gsl_expm1(x);
	}
	inline double gmath::hypot ( double x, double y )
	{
		return gsl_hypot(x,y);
	}
	inline double gmath::acosh ( double x )
	{
		return gsl_acosh(x);
	}
	inline double gmath::asinh ( double x )
	{
		return gsl_asinh(x);
	}
	inline double gmath::atanh ( double x )
	{
		return gsl_atanh(x);
	}
	inline double gmath::ldexp ( double x, int e )
	{
		return gsl_ldexp(x,e);
	}
	inline double gmath::frexp ( double x, int* e )
	{
		return gsl_frexp(x, e);
	}
	inline double gmath::pow_int ( double x, int n )
	{
		return gsl_pow_int(x,n);
	}
	inline double gmath::pow_2 ( double x )
	{
		return gsl_pow_2(x);
	}
	inline double gmath::pow_3 ( double x )
	{
		return gsl_pow_3(x);
	}
	inline double gmath::pow_4 ( double x )
	{
		return gsl_pow_4(x);
	}
	inline double gmath::pow_5 ( double x )
	{
		return gsl_pow_5(x);
	}
	inline double gmath::pow_6 ( double x )
	{
		return gsl_pow_6(x);
	}
	inline double gmath::pow_7 ( double x )
	{
		return gsl_pow_7(x);
	}
	inline double gmath::pow_8 ( double x )
	{
		return gsl_pow_8(x);
	}
	inline double gmath::pow_9 ( double x )
	{
		return gsl_pow_9(x);
	}
	inline int gmath::sign ( float x )
	{
		return GSL_SIGN(x);
	}
	inline int gmath::sign ( double x )
	{
		return GSL_SIGN(x);
	}
	inline int gmath::sign ( int x )
	{
		return GSL_SIGN(x);
	}
	inline int gmath::sign ( long x )
	{
		return GSL_SIGN(x);
	}
	inline bool gmath::is_odd ( int x )
	{
		return GSL_IS_ODD(x);
	}
	inline bool gmath::is_even ( int x )
	{
		return GSL_IS_EVEN(x);
	}
	inline int gmath::max ( int x, int y )
	{
		return GSL_MAX_INT(x,y);
	}
	inline double gmath::max ( double x, double y )
	{
		return GSL_MAX_DBL(x,y);
	}
	inline long double gmath::max ( long double x, long double y )
	{
		return GSL_MAX_LDBL(x,y);
	}
	inline int gmath::min ( int x, int y )
	{
		return GSL_MIN_INT(x,y);
	}
	inline double gmath::min ( double x, double y )
	{
		return GSL_MIN_DBL(x,y);
	}
	inline long double gmath::min ( long double x, long double y )
	{
		return GSL_MIN_LDBL(x,y);
	}
	inline int gmath::fcmp ( double x, double y, double epsilon )
	{
		return gsl_fcmp(x,y,epsilon);
	}
}

#endif   /* ----- #ifndef GMATH_INC  ----- */

